The second part of the program simulates the distribution of the t-statistic for the unit root test under both models. The results are first stored in a matrix object. Then the results are exported to a second workfile where the matrix results are converted to series for further processing. In particular, we plot the kernel density estimates of the simulated distributions of the t-statistic and compute the size and power of the ADF test.
' program to simulate the small sample distribution of ADF test
' for version 4
' last checked 11/2/2000 h
'
' Rudebusch (1993) "The Uncertain Unit Root in Real GNP," AER,
' 83, 1, 264-272.
' set up workfile
workfile adftest q 1948:1 1988:4
' fill GDP series
series gdp
gdp.fill 1284,1295.69995117,1303.80004883,1316.40002441,1305.30004883,1302,1312.59997559,1301.90002441,1350.90002441,1393.5,1445.19995117,1484.5,1504.09997559,1548.30004883,1585.40002441,1596,1607.69995117,1612.09997559,1621.90002441,1657.80004883,1687.30004883,1695.30004883,1687.90002441,1671.19995117,1660.80004883,1658.40002441,1677.69995117,1698.30004883,1742.5,1758.59997559,1778.19995117,1793.90002441,1787,1798.5,1802.19995117,1826.59997559,1836.40002441,1834.80004883,1851.19995117,1830.5,1790.09997559,1804.40002441,1840.90002441,1880.90002441,1904.90002441,1937.5,1930.80004883,1941.90002441,1976.90002441,1971.69995117,1973.69995117,1961.09997559,1977.40002441,2006,2035.19995117,2076.5,2103.80004883,2125.69995117,2142.60009766,2140.19995117,2170.89990234,2199.5,2237.60009766,2254.5,2311.10009766,2329.89990234,2357.39990234,2364,2410.10009766,2442.80004883,2485.5,2543.80004883,2596.80004883,2601.39990234,2626.10009766,2640.5,2657.19995117,2669,2699.5,2715.10009766,2752.10009766,2796.89990234,2816.80004883,2821.69995117,2864.60009766,2867.80004883,2884.5,2875.10009766,2867.80004883,2859.5,2895,2873.30004883,2939.89990234,2944.19995117,2962.30004883,2977.30004883,3037.30004883,3089.69995117,3125.80004883,3175.5,3253.30004883,3267.60009766,3264.30004883,3289.10009766,3259.39990234,3267.60009766,3239.10009766,3226.39990234,3154,3190.39990234,3249.89990234,3292.5,3356.69995117,3369.19995117,3381,3416.30004883,3466.39990234,3525,3574.39990234,3567.19995117,3591.80004883,3707,3735.60009766,3779.60009766,3780.80004883,3784.30004883,3807.5,3814.60009766,3830.80004883,3732.60009766,3733.5,3808.5,3860.5,3844.39990234,3864.5,3803.10009766,3756.10009766,3771.10009766,3754.39990234,3759.60009766,3783.5,3886.5,3944.39990234,4012.10009766,4089.5,4144,4166.39990234,4194.20019531,4221.79980469,4254.79980469,4309,4333.5,4390.5,4387.70019531,4412.60009766,4427.10009766,4460,4515.29980469,4559.29980469,4625.5,4655.29980469,4704.79980469,4734.5,4779.70019531
' work with log transformation
series ly = log(gdp)
' estimate trend stationary model
equation eq_ts.ls ly c @trend+1 ly(-1) ly(-2)
' estimate difference stationary model
equation eq_ds.ls d(ly) c d(ly(-1))
' compute impulse response from each model
' initialize impulse and response series
series u = 0
series ly_ts = 0
series ly_ds = 0
' create the impulse series U
smpl 48:3 48:3
series u = 1
smpl 48:3 88:4
' compute impulse response from trend stationary model
ly_ts = eq_ts.@coefs(3)*ly_ts(-1) + eq_ts.@coefs(4)*ly_ts(-2) + u
' compute impulse response from difference stationary model
ly_ds = ly_ds(-1) + eq_ds.@coefs(2)*d(ly_ds(-1)) + u
' plot the impulse responses
smpl 48:3 58:3
graph gh1.line ly_ts ly_ds
' label graph
gh1.setelem(1) legend(Trend stationary model)
gh1.setelem(2) legend(Difference stationary model)
gh1.addtext(t) Impulse Response Functions
show gh1
'--------------------------------------------------------------------------
' simulate ADF t-stats from each model
'--------------------------------------------------------------------------
' initialize simulated series
smpl 48:1 48:2
series lyts = ly
series lyds = ly
' declare matrix to store simulated ADF t-stats
!n = 1000 ' number of replications
matrix(!n,2) adft
' declare working equation
equation eq_test
' monte carlo loop
rndseed 123456789
smpl 48:3 @last
for !i = 1 to !n
' trend stationary model
' simulate series
lyts = eq_ts.@coefs(1) + eq_ts.@coefs(2)*(@trend+1) + eq_ts.@coefs(3)*lyts(-1) + eq_ts.@coefs(4)*lyts(-2) + eq_ts.@se*nrnd
' run ADF test
eq_test.ls d(lyts) lyts(-1) c @trend+1 d(lyts(-1))
' store t-stat
adft(!i,1) = eq_test.@tstat(1)
' difference stationary model
' simulate series
lyds = lyds(-1) + eq_ds.@coefs(1) + eq_ds.@coefs(2)*d(lyds(-1)) + eq_ds.@se*nrnd
' run ADF test
eq_test.ls d(lyds) lyds(-1) c @trend+1 d(lyds(-1))
' store t-stat
adft(!i,2) = eq_test.@tstat(1)
next
' store actual ADF t-stat
eq_test.ls d(ly) ly(-1) c @trend+1 d(ly(-1))
!tadf = eq_test.@tstat(1)
' store results to database
db db_mcarlo
store adft
' close workfile
close adftest
'--------------------------------------------------------------------------
' process results as series in new workfile
'--------------------------------------------------------------------------
' create new workfile
workfile result u 1 !n
' fetch results
fetch adft
' declare series
series t_ts
series t_ds
' convert simulated t-stats to series
group gr1 t_ts t_ds
mtos(adft,gr1)
' determine knots for kernel density estimates
' to ensure it fits in current workfile
!knots = 100
if (!knots>!n) then
!knots = !n
endif
smpl 1 !knots
for %mod ts ds
' store kernel density estimates
do t_{%mod}.kdensity(!knots, o=kdmat_{%mod})
' convert density estimates to series
series kd{%mod}_x
series kd{%mod}_f
group g{%mod} kd{%mod}_x kd{%mod}_f
mtos(kdmat_{%mod}, g{%mod})
next
'compute empirical size of test
' area of H0 to the left of !tadf
smpl @all if t_ds<=!tadf
!siz = @obs(t_ds)/!n
'compute empirical power of test
' area of H1 to the left of !tadf
smpl @all if t_ts<=!tadf
!pow = @obs(t_ts)/!n
' plot the two density estimates in one graph
smpl 1 !knots
group g3 gts gds
freeze(graph1) g3.xyline(b)
'graph1.setelem(1) legend(Alternative (trend stationary))
'graph1.setelem(2) legend(Null (difference stationary))
graph1.addtext(t) Kernel density estimates of simulated t-distribution
graph1.legend -display
graph1.addtext(1.2,0.1) H1 (trend stationary)
graph1.addtext(1.8,0.8) H0 (difference stationary)
' draw vertical dashline at actual t-stat
graph1.draw(dashline, bottom, rgb(156,156,156)) !tadf
graph1.addtext(2.8,2.4) size = !siz
graph1.addtext(2.8,2.6) power = !pow
show graph1